Species data cleaning
Florencia Grattarola
30/04/2022
library(patchwork)
library(rnaturalearth)
library(countrycode)
library(lubridate)
library(tmap)
library(leaflet)
library(terra)
library(sf)
sf::sf_use_s2(FALSE)
library(tidyverse)1 Basemaps
# world and Latin America
world <- rnaturalearth::ne_countries(scale = 'large', returnclass = 'sf')
bbox_Latam_unprojected <- c(xmin=-118.40137, ymin=-55.89170, xmax=-34.80547, ymax= 32.71533)
Latam_unprojected <- world %>% st_crop(bbox_Latam_unprojected)
# equal area projection (Equatorial Lambert azimuthal equal-area)
equalareaCRS <- '+proj=laea +lon_0=-73.125 +lat_0=0 +datum=WGS84 +units=m +no_defs'
Latam_projected <-st_read('data/Latam_vector.shp')## Reading layer `Latam_vector' from data source
## `/Users/flograttarola/Documents/GitHub/yaguarundi_IDM/data/Latam_vector.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -4408686 ymin: -6083095 xmax: 4161425 ymax: 3827660
## CRS: unknownLatam_countries <- sf::st_read('data/Latam_vector_countries.shp', quiet = T) %>% sf::st_transform(crs=equalareaCRS)
Latam <- st_union(st_make_valid(Latam_projected))
# raster of Latin America to use as template
# dimension of raster layers 100km^2 (100000m^2)
Latam.raster <- terra::rast('data/Latam_raster.tif')
# species list
# SOURCE: Mammal Diversty Database AMS
speciesList_MDD <- read_csv('data/MDD_Carnivora_List.csv')
# example species selected
selectedSpecies <- c('Leopardus geoffroyi',
'Herpailurus yagouaroundi',
'Lontra longicaudis',
'Tremarctos ornatus',
'Cerdocyon thous',
'Eira barbara')# IUCN range distribution
yaguarundi_IUCN <- sf::st_read('data/yaguarundi_IUCN.shp', quiet = T) %>% sf::st_transform(crs=equalareaCRS)
MAP.IUCN <- tm_graticules(alpha = 0.3) +
tm_shape(Latam_countries) +
tm_fill(col='grey90') +
tm_borders(col='grey60', alpha = 0.4) +
tm_shape(yaguarundi_IUCN) +
tm_fill(col='grey20', alpha = 0.4) +
tm_layout(legend.outside = T, frame.lwd = 0.3, scale=1.2, legend.outside.size = 0.1)
MAP.IUCN
# tmap_save(MAP.IUCN, filename='docs/figs/MAP.IUCN.svg', dpi=300, height = 7, width = 6, device = svglite::svglite)2 Data
2.1 Presence/absence and abundance data
Terms standardised:
- recordID
- species
- temporalSpan
- yearEnd
- decimalLatitude
- decimalLongitude
- area_ha (study area)
- area_m2
- countryCode
- stateProvince
- locality
- presence
- abundance
- effort
- effortUnit
- dataset
pointsAndSurveys_NeotropCarnivores <-
read_csv('data/NEOTROPICAL_CARNIVORES_DATASET_2020-04.csv',
guess_max = 90000) # Neotropical carnivores Data Paper
data_PA_Ab <-
pointsAndSurveys_NeotropCarnivores %>%
# filter the records to include only the Neotropical Mammals
filter(SPECIES %in% speciesList_MDD$sp) %>%
# filter only the records taken with camera traps -> assumption: these will reflect presence-absece data
filter(METHOD=='Camera trap') %>%
# filter records with no coordinates
filter(!is.na(LAT_Y)) %>%
# filter occurrences that don't have month and year when they were recorded
filter(!is.na(COL_END_YR)&!is.na(COL_START_MO)&!is.na(COL_END_MO)&!is.na(COL_START_YR)) %>%
# one area has a range instead of a value, I took the mean of the range as area
mutate(AREA_HA=ifelse(AREA_HA=='700-1000', mean(c(700, 1000)), AREA_HA)) %>%
# filter records with no area of study recorded / or with sampling level area and precision recorded
filter(!is.na(AREA_HA) | (is.na(AREA_HA) & !is.na(PRECISION_m) & SAMPLING_LEVEL=='AREA')) %>%
mutate(area_ha=as.numeric(AREA_HA),
area_m2=ifelse(!is.na(AREA_HA),
area_ha*10000, # 1 hectare are 10000 meters / will be used for buffers
PRECISION_m)) %>%
# to have the same values as the presence-only I kept country codes
mutate(countryCode=countrycode::countrycode(COUNTRY,
origin = 'country.name',
destination = 'iso2c')) %>%
# since the dataset has no day, dates were assumed to start and end the first day of the month
mutate(date_start=lubridate::dmy(str_c('1', COL_START_MO, COL_START_YR))) %>%
mutate(date_end=lubridate::dmy(str_c('1', COL_END_MO, COL_END_YR))) %>%
# the temporal span was calculated -> there are errors in some dates, and the temporal span is negative
mutate(temporalSpan=ifelse(date_end-date_start<0, date_start-date_end, date_end-date_start)) %>%
# decimal places for some records corrected
mutate(EFFORT=ifelse(grepl('\\.', EFFORT) & (temporalSpan>=31 | temporalSpan>=30),
str_remove(EFFORT, '\\.' ), EFFORT)) %>%
# the sampling effort was standardized to days
mutate(samplingEffort=ifelse(grepl('hours', EFFORT_UNIT, ignore.case = T), as.numeric(EFFORT)/24,
ifelse(grepl('days', EFFORT_UNIT, ignore.case = T), as.numeric(EFFORT), NA))) %>%
# filter those without sampling effort (no data was provided and it cannot be derived)
filter(!is.na(samplingEffort)) %>%
mutate(effortUnit='days') %>%
# calculate the number of camera traps needed
mutate(numCameras=floor(samplingEffort/temporalSpan)) %>%
# filter unreliable values
mutate(temporalSpan=ifelse(numCameras==Inf & temporalSpan==0 & samplingEffort>31, NA,
ifelse(numCameras==Inf & temporalSpan==0 & samplingEffort<=31,
floor(as.numeric(EFFORT)), temporalSpan))) %>%
filter(!is.na(temporalSpan)) %>%
mutate(independentLocation=str_c(LAT_Y,':',LONG_X)) %>%
mutate(independentYearSpan=str_c(COL_START_YR,':',COL_END_YR)) %>%
#mutate(yearEnd=lubridate::year(date_end)) %>%
dplyr::select(recordID=ID,
date_start, date_end,
species=SPECIES,
temporalSpan,
#yearEnd,
independentLocation,
independentYearSpan,
decimalLatitude=LAT_Y,
decimalLongitude=LONG_X,
area_ha,
area_m2,
countryCode,
stateProvince=STATE,
locality=SITE,
presence=OCCUR,
abundance=N_RECORDS, # not all of the records have abundance data
effort=samplingEffort,
effortUnit,
dataset=DATASET)%>%
mutate(dataset='pointsAndSurveys_NeotropCarnivores') # in total 34,389 records
# Generate 0 (zeros) for locations where the species hasn't been recorded
# this can very probably be done in a more efficient way, but this one works!
data_PA_Ab_0 <- data_PA_Ab %>%
pivot_wider(names_from = species,
values_from = c(presence, abundance, effort),
values_fill = c(list(presence=0),
list(abundance=0),
list(effort=0))) %>%
pivot_longer(cols=starts_with(c('presence_', 'abundance_', 'effort_')),
names_to=c('metric', 'species'),
names_sep='_') %>%
pivot_wider(names_from = metric,
values_from = value,
values_fill = c(list(value=0))) %>%
distinct(species, independentYearSpan, independentLocation, presence, .keep_all = T) %>%
group_by(species, independentYearSpan, independentLocation) %>%
mutate(presence=sum(presence, na.rm = T),
abundance=sum(abundance, na.rm = T),
effort=max(effort, na.rm = T)) %>%
distinct(species, independentYearSpan, independentLocation, .keep_all = T) %>%
ungroup() %>% group_by(independentLocation, independentYearSpan) %>%
mutate(effort=max(effort, na.rm = T)) %>%
ungroup()2.2 Presence-only data
Terms standardised:
- recordID,
- species,
- year,
- eventDate,
- decimalLatitude,
- decimalLongitude,
- countryCode, - stateProvince,
- locality,
- dataset
# The data download and cleaning can be found in the GBIF_download.R script
occurrencePoints_GBIF <- read_csv('data/0046613-210914110416597_CLEAN.csv',
guess_max = 120000) # GBIF data (cleaned)
data_PO <-
occurrencePoints_GBIF %>%
# rename the species yaguarundi for GBIF data
mutate(species=ifelse(species=='Puma yagouaroundi', 'Herpailurus yagouaroundi', species)) %>%
# filter records from unrealistic or unwanted) locations: Japan, USA and Canada
filter(countryCode!='JP' & countryCode!='US' & countryCode!='CA') %>%
# filter occurrences that don't have the year when they were recorded
filter(!is.na(year) & year<=2020) %>%
# transform eventdate to class datetime
mutate(eventDate=lubridate::ymd_hms(eventDate)) %>%
# remove duplicates according to species, date, latitude and longitude
group_by(species, eventDate, decimalLatitude, decimalLongitude) %>%
slice_head(n=1) %>% ungroup() %>%
# select columns of interest
dplyr::select(recordID=gbifID,
species,
year, eventDate,
decimalLatitude,
decimalLongitude,
countryCode,
stateProvince,
locality) %>%
mutate(dataset='occurrencePoints_GBIF') 2.3 Check data for each years on example species
# Presence-only
data_PO %>% filter(species == 'Herpailurus yagouaroundi') %>% filter(year>=2000) %>%
mutate(period=ifelse(year<=2014, 'pre2014', 'post2014')) %>%
group_by(species, period) %>% count()## # A tibble: 2 × 3
## # Groups: species, period [2]
## species period n
## <chr> <chr> <int>
## 1 Herpailurus yagouaroundi post2014 234
## 2 Herpailurus yagouaroundi pre2014 193# Presence-absence
data_PA_Ab_0 %>% filter(species %in% selectedSpecies & presence ==1) %>%
mutate(period=ifelse(year(date_start)<=2014, 'pre2014', 'post2014')) %>%
group_by(species, period) %>% count()## # A tibble: 12 × 3
## # Groups: species, period [12]
## species period n
## <chr> <chr> <int>
## 1 Cerdocyon thous post2014 874
## 2 Cerdocyon thous pre2014 1232
## 3 Eira barbara post2014 595
## 4 Eira barbara pre2014 1198
## 5 Herpailurus yagouaroundi post2014 215
## 6 Herpailurus yagouaroundi pre2014 399
## 7 Leopardus geoffroyi post2014 70
## 8 Leopardus geoffroyi pre2014 59
## 9 Lontra longicaudis post2014 41
## 10 Lontra longicaudis pre2014 30
## 11 Tremarctos ornatus post2014 17
## 12 Tremarctos ornatus pre2014 63 Transform data for the model integration
3.1 Presence-absence data
To assess the presences and absences we created a buffer area around the study location according to the area reported. Then we combined the overlapping buffers into blobs.
Notes from Nagy-Reis et al. (2020): For data labeled as SAMPLING_LEVEL= AREA, the PRECISION_m field could be used to get a sense of the study area: Page 45: “The centroid data are identified by having”AREA” in the ‘SAMPLING_LEVEL’ attribute, which indicates that the exact location is unknown and its precision reflects the size of study area (‘AREA_HA’).” But if SAMPLING_LEVEL = UNIT, then PRECISION_m no longer informs on area size, it simply informs on the precision of the method used to collect the coordinate and I would personally not recommend using PRECISION_m in that case.
3.1.1 Create blobs.
Blobs contain data on total effort, total temporal span, a list of the records IDs that are included for the blob and an area size of the combined areas.
data_PA_Ab_0_GIS <- data_PA_Ab_0 %>%
filter(species == 'Herpailurus yagouaroundi') %>%
mutate(year=lubridate::year(date_end)) %>%
filter(year>=2000, year<=2021) %>%
as.data.frame() %>%
sf::st_as_sf(coords=c('decimalLongitude', 'decimalLatitude')) %>%
st_set_crs(4326)
PA <- st_transform(data_PA_Ab_0_GIS, crs=equalareaCRS)
# create a blob for all sites
# blobs are definded as areas where some sampling was done over the entire period
periodStart <- ymd('2000-01-01')
periodEnd <- ymd('2014-01-01')
periodMax <- ymd('2020-12-31')
PA_allsites <- PA %>%
mutate(start_date=as_date(ifelse(date_end-date_start<0, date_end, date_start)),
end_date=as_date(ifelse(date_end-date_start<0, date_start, date_end))) %>%
mutate(span=interval(date_start, date_end),
target= interval(periodEnd, today()),
period=ifelse(span %within% target, 'pos', 'pre')) %>%
distinct(independentLocation, period, .keep_all = T)
PA_allsites_buff <- st_buffer(PA_allsites, sqrt(PA_allsites$area_m2/pi))
blobs_pre <- PA_allsites_buff %>% filter(period=='pre') %>%
st_union() %>% st_cast('POLYGON') %>% st_sf('ID'= 1:length(.))
blobs_pos <- PA_allsites_buff %>% filter(period=='pos') %>%
st_union() %>% st_cast('POLYGON') %>% st_sf('ID'= 1:length(.))
# to calculate the effort in days for each blob, we need to differenciate the two periods
PA_allsites_effort_and_time <- PA %>%
mutate(start_date=as_date(ifelse(date_end-date_start<0, date_end, date_start)),
end_date=as_date(ifelse(date_end-date_start<0, date_start, date_end))) %>%
mutate(span=interval(date_start, date_end),
target= interval(periodEnd, today()),
period=ifelse(span %within% target, 'pos', 'pre')) %>%
group_by(independentLocation, independentYearSpan, area_m2) %>%
summarise(effort=max(effort),
period=paste(unique(period)),
recordIDs=paste(unique(recordID), collapse = ';'),
maxEndDate=max(end_date),
minStartDate=min(start_date)) %>%
mutate(maxTemporalSpan= maxEndDate-minStartDate)
# to calculate the effort in days for each blob, we need to differenciate the two periods
blobs_efforttime_pre <- st_join(blobs_pre,
PA_allsites_effort_and_time %>% filter(period=='pre'),
left=TRUE) %>%
group_by(ID) %>%
summarise(effort=sum(effort),
temporalSpan=sum(maxTemporalSpan),
recordIDs=paste(unique(recordIDs), collapse = ';')) %>%
mutate(blobArea=st_area(.))
blobs_efforttime_pos <- st_join(blobs_pos,
PA_allsites_effort_and_time %>% filter(period=='pos'),
left=TRUE) %>%
group_by(ID) %>%
summarise(effort=sum(effort),
temporalSpan=sum(maxTemporalSpan),
recordIDs=paste(unique(recordIDs), collapse = ';')) %>%
mutate(blobArea=st_area(.))The data for the selected species
# function to create blobs for pre and pos 2010 for each species
calculate_PA_sp_blobs <- function(PA, blobs, sp, prepos){
periodStart <- ymd('2000-01-01')
periodEnd <- ymd('2014-01-01')
df_period <- PA %>% filter(species==sp & presence==1) %>%
filter(date_start>=periodStart) %>%
mutate(span=interval(date_start, date_end),
target= interval(periodEnd, today()),
period=ifelse(span %within% target, 'pos', 'pre')) %>%
filter(period==prepos) %>% select(species, presence)
df_period_blobs <- st_join(blobs, df_period,
left=TRUE, join = st_contains) %>%
group_by(ID) %>%
summarise(presence=max(presence),
temporalSpan=max(temporalSpan),
effort=max(effort),
blobArea=max(blobArea)) %>%
mutate(presence=ifelse(is.na(presence), 0, presence))
return(df_period_blobs)
}
# Herpailurus yagouaroundi
PA_herpailurus_pre_blobs <- calculate_PA_sp_blobs(PA, blobs_efforttime_pre, 'Herpailurus yagouaroundi', 'pre')
PA_herpailurus_pos_blobs <- calculate_PA_sp_blobs(PA, blobs_efforttime_pos, 'Herpailurus yagouaroundi', 'pos')3.1.2 Quick check
# quick look un how data is summarised
st_join(blobs_efforttime_pos, PA %>% filter(species=='Herpailurus yagouaroundi' & presence==1) %>%
mutate(span=interval(date_start, date_end),
target= interval(ymd('2014-01-01'), today()),
period=ifelse(span %within% target, 'pos', 'pre')) %>%
filter(period=='pos') %>% select(species, presence),
left=TRUE, join = st_contains) %>%
group_by(ID) %>%
summarise(presence=max(presence),
temporalSpan=max(temporalSpan),
effort=max(effort),
blobArea=max(blobArea)) %>%
mutate(presence=ifelse(is.na(presence), 0, presence)) %>% head()## Simple feature collection with 6 features and 5 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 2146306 ymin: -3500274 xmax: 2316801 ymax: -3237134
## CRS: +proj=laea +lon_0=-73.125 +lat_0=0 +datum=WGS84 +units=m +no_defs
## # A tibble: 6 × 6
## ID presence temporalSpan effort blobArea .
## <int> <dbl> <drtn> <dbl> [m^2] <POLYGON [m]>
## 1 1 0 224 days 2464 2398303. ((2147016 -3500112, 2147004 -35…
## 2 2 0 372 days 3752 2992037. ((2178918 -3462334, 2178899 -34…
## 3 3 0 960 days 6592 89167699. ((2179573 -3381373, 2179413 -33…
## 4 4 0 1782 days 18372 2636198. ((2159621 -3324197, 2159597 -33…
## 5 5 1 457 days 1191 15992690. ((2264845 -3282291, 2264773 -32…
## 6 6 1 670 days 1640 468832531. ((2304059 -3258568, 2303746 -32…3.2 Presence-only data
# rasterisation of the presence-only data
data_PO_GIS <- data_PO %>%
filter(species == 'Herpailurus yagouaroundi') %>%
mutate(year=lubridate::year(eventDate)) %>%
filter(year>=2000, year<=2021) %>%
as.data.frame() %>%
sf::st_as_sf(coords=c('decimalLongitude', 'decimalLatitude')) %>%
st_set_crs(4326)
PO <- st_transform(data_PO_GIS, crs=equalareaCRS)
# functions to create presence-only rasters for pre and pos 2014 fo each species
calculate_PO_sp_period_raster <- function(PO, sp, raster, prepos){
periodStart <- 2000
periodEnd <- 2014
df_period <- PO %>% filter(species==sp & !is.na(year)) %>%
filter(year>=periodStart) %>%
mutate(period=ifelse(year>=periodEnd , 'pos', 'pre'),
occurrence=1) %>% filter(period==prepos)
df_period_raster <- terra::rasterize(x = vect(df_period),
y = raster,
field = 'occurrence',
fun = 'length',
sum = T, background=0) %>% mask(., raster)
df_country_raster <- terra::rasterize(x = vect(df_period),
y = raster,
field = 'countryCode',
fun = 'first') %>% mask(., raster)
names(df_period_raster) <- 'count'
names(df_country_raster) <- 'country'
df_raster <- c(df_period_raster, df_country_raster)
return(df_raster)
}
# Herpailurus yagouaroundi
PO_herpailurus_pre_raster <- calculate_PO_sp_period_raster(PO, 'Herpailurus yagouaroundi', Latam.raster, 'pre')
PO_herpailurus_pos_raster <- calculate_PO_sp_period_raster(PO, 'Herpailurus yagouaroundi', Latam.raster, 'pos')3.3 Data visualisation
3.3.1 Manuscript plots
# plot presence-absence data
data_PA_pre_pos <- rbind(PA_herpailurus_pre_blobs %>% mutate(period='pre'),
PA_herpailurus_pos_blobs %>% mutate(period='pos'))
gg_PA_pre <- ggplot() +
geom_sf(data=Latam_countries, fill='#fafaf8', col='#bfc5c7', size=0.2) +
geom_sf(data=yaguarundi_IUCN, fill='grey60', alpha=0.5, col='grey90', size=0.1) +
geom_sf(data=data_PA_pre_pos %>% filter(period=='pre'),
aes(fill=factor(presence), col=factor(presence)), size=1, show.legend = F) +
coord_sf(xlim = c(-3500000, 4100000), ylim = c(-4400000, 3000000)) + # st_bbox(yaguarundi_IUCN)
scale_fill_manual(values = c("#474545","#E31A1C"))+
scale_color_manual(values = c("#474545","#E31A1C")) +
labs(title='PA', subtitle='pre: 2000-2013', col='') +
theme_bw() +
theme(text=element_text(size = 14))
gg_PA_pos <- ggplot() +
geom_sf(data=Latam_countries, fill='#fafaf8', col='#bfc5c7', size=0.2) +
geom_sf(data=yaguarundi_IUCN, fill='grey60', alpha=0.5, col='grey90', size=0.1) +
geom_sf(data=data_PA_pre_pos %>% filter(period=='pos'),
aes(fill=factor(presence), col=factor(presence)), size=1, show.legend = F) +
coord_sf(xlim = c(-3500000, 4100000), ylim = c(-4400000, 3000000)) + # st_bbox(yaguarundi_IUCN)
scale_fill_manual(values = c("#434445","#1F78B4"))+
scale_color_manual(values = c("#434445","#1F78B4")) +
labs(title='', subtitle='pos: 2014 to 2021', col='') +
theme_bw() +
theme(text=element_text(size = 14))
gg_PA_pre + gg_PA_pos
# gg_PA <- gridExtra::grid.arrange(gg_PA_pre, gg_PA_pos, nrow=1)
# ggsave(plot = gg_PA, filename='docs/figs/presence-absence-plots.svg', device = 'svg', width=10, height=7, dpi=300)data_PO_pre_pos <- data_PO_GIS %>%
mutate(period=ifelse(year>=2014 , 'pos', 'pre'), occurrence=1)
gg_PO_pre <- ggplot() +
geom_sf(data=Latam_countries, fill='#fafaf8', col='#bfc5c7', size=0.2) +
geom_sf(data=yaguarundi_IUCN, fill='grey60', alpha=0.5, col='grey90', size=0.1) +
geom_sf(data=data_PO_pre_pos %>% filter(period=='pre'), col="#E31A1C", size=0.5, show.legend = F) +
coord_sf(xlim = c(-3500000, 4100000), ylim = c(-4400000, 3000000)) + # st_bbox(yaguarundi_IUCN)
labs(title='PO pre: 2000-2013', col='') +
theme_bw() +
theme(text=element_text(size = 14))
gg_PO_pos <- ggplot() +
geom_sf(data=Latam_countries, fill='#fafaf8', col='#bfc5c7', size=0.2) +
geom_sf(data=yaguarundi_IUCN, fill='grey60', alpha=0.5, col='grey90', size=0.1) +
geom_sf(data=data_PO_pre_pos %>% filter(period=='pos'), col="#1F78B4", size=0.5, show.legend = F) +
coord_sf(xlim = c(-3500000, 4100000), ylim = c(-4400000, 3000000)) + # st_bbox(yaguarundi_IUCN)
labs(title='PO post: 2014 to 2021', col='') +
theme_bw() +
theme(text=element_text(size = 14))
gg_PO_pre + gg_PO_pos
# gg_PO <- gridExtra::grid.arrange(gg_PO_pre, gg_PO_pos, nrow=1)
# ggsave(plot = gg_PO, filename='docs/figs/presence-only-plots.svg', device = 'svg', width=10, height=7, dpi=300)3.3.2 Interactive plot
Species Herpailurus yagouaroundi, POST period (2014-2020). Blobs of presence/absence and presence-only records.
In red the blobs (polygons) where the species is absent, in white where the species is present. Blue dots are the lat/long locations of each study where the species was present. Black dots are presence-only records.
4 Save data
# blobs
saveRDS(blobs_efforttime_pre, 'data/blobs_pre.Rds')
saveRDS(blobs_efforttime_pos, 'data/blobs_pos.Rds')
# presence-only rasters
terra::writeRaster(PO_herpailurus_pre_raster, 'data/PO_herpailurus_pre_raster.tif', overwrite=T)
terra::writeRaster(PO_herpailurus_pos_raster, 'data/PO_herpailurus_pos_raster.tif', overwrite=T)
# presences/absences blobs / polygons
saveRDS(PA_herpailurus_pre_blobs, 'data/PA_herpailurus_pre_blobs.Rds')
saveRDS(PA_herpailurus_pos_blobs, 'data/PA_herpailurus_pos_blobs.Rds')